function [c] = polynomial(x,y,z);
%POLYNOMIAL  Fit a 5th order polynomial in 3D to a given set of points

% Constant is assumed to be 1
x0y0z0 = ones(56,1);

% Compute all the polynomial values at each point
x1y0z0 = x;
x2y0z0 = x1y0z0 .* x1y0z0;
x3y0z0 = x2y0z0 .* x1y0z0;
x4y0z0 = x3y0z0 .* x1y0z0;
x5y0z0 = x4y0z0 .* x1y0z0;

x0y1z0 = y;
x1y1z0 = x1y0z0 .* x0y1z0;
x2y1z0 = x2y0z0 .* x0y1z0;
x3y1z0 = x3y0z0 .* x0y1z0;
x4y1z0 = x4y0z0 .* x0y1z0;

x0y2z0 = x0y1z0 .* x0y1z0;
x1y2z0 = x1y1z0 .* x0y1z0;
x2y2z0 = x2y1z0 .* x0y1z0;
x3y2z0 = x3y1z0 .* x0y1z0;

x0y3z0 = x0y2z0 .* x0y1z0;
x1y3z0 = x1y2z0 .* x0y1z0;
x2y3z0 = x2y2z0 .* x0y1z0;

x0y4z0 = x0y3z0 .* x0y1z0;
x1y4z0 = x1y3z0 .* x0y1z0;

x0y5z0 = x0y4z0 .* x0y1z0;

x0y0z1 = z;
x1y0z1 = x0y0z1 .* x1y0z0;
x2y0z1 = x1y0z1 .* x1y0z0;
x3y0z1 = x2y0z1 .* x1y0z0;
x4y0z1 = x3y0z1 .* x1y0z0;

x0y1z1 = x0y1z0 .* x0y0z1;
x1y1z1 = x1y1z0 .* x0y0z1;
x2y1z1 = x2y1z0 .* x0y0z1;
x3y1z1 = x3y1z0 .* x0y0z1;

x0y2z1 = x0y2z0 .* x0y0z1;
x1y2z1 = x1y2z0 .* x0y0z1;
x2y2z1 = x2y2z0 .* x0y0z1;

x0y3z1 = x0y3z0 .* x0y0z1;
x1y3z1 = x1y3z0 .* x0y0z1;

x0y4z1 = x0y4z0 .* x0y0z1;

x0y0z2 = x0y0z1 .* x0y0z1;
x1y0z2 = x1y0z1 .* x0y0z1;
x2y0z2 = x2y0z1 .* x0y0z1;
x3y0z2 = x3y0z1 .* x0y0z1;

x0y1z2 = x0y1z1 .* x0y0z1;
x1y1z2 = x1y1z1 .* x0y0z1;
x2y1z2 = x2y1z1 .* x0y0z1;

x0y2z2 = x0y2z1 .* x0y0z1;
x1y2z2 = x1y2z1 .* x0y0z1;

x0y3z2 = x0y3z1 .* x0y0z1;

x0y0z3 = x0y0z2 .* x0y0z1;
x1y0z3 = x1y0z2 .* x0y0z1;
x2y0z3 = x2y0z2 .* x0y0z1;

x0y1z3 = x0y1z2 .* x0y0z1;
x1y1z3 = x1y1z2 .* x0y0z1;

x0y2z3 = x0y2z2 .* x0y0z1;

x0y0z4 = x0y0z3 .* x0y0z1;
x1y0z4 = x1y0z3 .* x0y0z1;

x0y1z4 = x0y1z3 .* x0y0z1;

x0y0z5 = x0y0z4 .* x0y0z1;

% Set up the linear system to solve for the polynomial coefficients
A = [x1y0z0, ...
     x2y0z0, ...
     x3y0z0, ...
     x4y0z0, ...
     x5y0z0, ...
     ...
     x0y1z0, ...
     x1y1z0, ...
     x2y1z0, ...
     x3y1z0, ...
     x4y1z0, ...
     ...
     x0y2z0, ...
     x1y2z0, ...
     x2y2z0, ...
     x3y2z0, ...
     ...
     x0y3z0, ...
     x1y3z0, ...
     x2y3z0, ...
     ...
     x0y4z0, ...
     x1y4z0, ...
     ...
     x0y5z0, ...
     ...
     x0y0z1, ...
     x1y0z1, ...
     x2y0z1, ...
     x3y0z1, ...
     x4y0z1, ...
     ...
     x0y1z1, ...
     x1y1z1, ...
     x2y1z1, ...
     x3y1z1, ...
     ...
     x0y2z1, ...
     x1y2z1, ...
     x2y2z1, ...
     ...
     x0y3z1, ...
     x1y3z1, ...
     ...
     x0y4z1, ...
     ...
     x0y0z2, ...
     x1y0z2, ...
     x2y0z2, ...
     x3y0z2, ...
     ...
     x0y1z2, ...
     x1y1z2, ...
     x2y1z2, ...
     ...
     x0y2z2, ...
     x1y2z2, ...
     ...
     x0y3z2, ...
     ...
     x0y0z3, ...
     x1y0z3, ...
     x2y0z3, ...
     ...
     x0y1z3, ...
     x1y1z3, ...
     ...
     x0y2z3, ...
     ...
     x0y0z4, ...
     x1y0z4, ...
     ...
     x0y1z4, ...
     ...
     x0y0z5  ...
     ];

% Solve the linear system to solve for the polynomial coefficients
c = A \ (-x0y0z0);
